CSCI E-25¶

Optical Flow¶

Steve Elston¶

In this assignment you will perform a exercises concerning optical flow. Code is provided and you are to answer the questions in the exercises.

These exercises are derived from an example in the Scikit Image User Guide.

Note: Before you proceed you may need to update the version of scikit-image and the Pooch package. To do so, uncomment and execute the code in the cell below. Then shut down your Jupyter server and restart.

In [1]:
!pip install --upgrade scikit-image
!pip install pooch

Requirement already satisfied: scikit-image in /opt/anaconda3/lib/python3.11/site-packages (0.22.0)
Collecting scikit-image
  Downloading scikit_image-0.23.2-cp311-cp311-macosx_12_0_arm64.whl.metadata (14 kB)
Requirement already satisfied: numpy>=1.23 in /opt/anaconda3/lib/python3.11/site-packages (from scikit-image) (1.26.4)
Requirement already satisfied: scipy>=1.9 in /opt/anaconda3/lib/python3.11/site-packages (from scikit-image) (1.11.4)
Requirement already satisfied: networkx>=2.8 in /opt/anaconda3/lib/python3.11/site-packages (from scikit-image) (3.1)
Requirement already satisfied: pillow>=9.1 in /opt/anaconda3/lib/python3.11/site-packages (from scikit-image) (10.2.0)
Requirement already satisfied: imageio>=2.33 in /opt/anaconda3/lib/python3.11/site-packages (from scikit-image) (2.33.1)
Requirement already satisfied: tifffile>=2022.8.12 in /opt/anaconda3/lib/python3.11/site-packages (from scikit-image) (2023.4.12)
Requirement already satisfied: packaging>=21 in /opt/anaconda3/lib/python3.11/site-packages (from scikit-image) (23.1)
Collecting lazy-loader>=0.4 (from scikit-image)
  Downloading lazy_loader-0.4-py3-none-any.whl.metadata (7.6 kB)
Downloading scikit_image-0.23.2-cp311-cp311-macosx_12_0_arm64.whl (13.3 MB)
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 13.3/13.3 MB 41.0 MB/s eta 0:00:0000:0100:01
Downloading lazy_loader-0.4-py3-none-any.whl (12 kB)
Installing collected packages: lazy-loader, scikit-image
  Attempting uninstall: lazy-loader
    Found existing installation: lazy_loader 0.3
    Uninstalling lazy_loader-0.3:
      Successfully uninstalled lazy_loader-0.3
  Attempting uninstall: scikit-image
    Found existing installation: scikit-image 0.22.0
    Uninstalling scikit-image-0.22.0:
      Successfully uninstalled scikit-image-0.22.0
Successfully installed lazy-loader-0.4 scikit-image-0.23.2
Collecting pooch
  Downloading pooch-1.8.1-py3-none-any.whl.metadata (9.5 kB)
Requirement already satisfied: platformdirs>=2.5.0 in /opt/anaconda3/lib/python3.11/site-packages (from pooch) (3.10.0)
Requirement already satisfied: packaging>=20.0 in /opt/anaconda3/lib/python3.11/site-packages (from pooch) (23.1)
Requirement already satisfied: requests>=2.19.0 in /opt/anaconda3/lib/python3.11/site-packages (from pooch) (2.31.0)
Requirement already satisfied: charset-normalizer<4,>=2 in /opt/anaconda3/lib/python3.11/site-packages (from requests>=2.19.0->pooch) (2.0.4)
Requirement already satisfied: idna<4,>=2.5 in /opt/anaconda3/lib/python3.11/site-packages (from requests>=2.19.0->pooch) (3.4)
Requirement already satisfied: urllib3<3,>=1.21.1 in /opt/anaconda3/lib/python3.11/site-packages (from requests>=2.19.0->pooch) (2.0.7)
Requirement already satisfied: certifi>=2017.4.17 in /opt/anaconda3/lib/python3.11/site-packages (from requests>=2.19.0->pooch) (2024.2.2)
Downloading pooch-1.8.1-py3-none-any.whl (62 kB)
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 63.0/63.0 kB 2.1 MB/s eta 0:00:00
Installing collected packages: pooch
Successfully installed pooch-1.8.1

To get started, execute the code in the cell below to import the packages you will need for this notebook.

In [1]:
import numpy as np
from matplotlib import pyplot as plt
from skimage.color import rgb2gray
from skimage.data import stereo_motorcycle, vortex
from skimage.transform import warp
from skimage.registration import optical_flow_tvl1, optical_flow_ilk

Iterative Lucas-Kanade (iLK) Algorithm¶

As a first simple example, you will apply the iterative Lucas-Kanade algorithm to a random field image. The random field image is generated by computing random pixel values. A second image is created by rotating the first around an axis in the image field. The goal is to register these two images and compute and display the optical flow field.

Now, execute the code in the cell below to load and display the two images.

In [2]:
def plot_grayscale(img, h=8, ax=None, title = ''):
    if(ax==None): fig, ax = plt.subplots(1,1, figsize=(h, h))
    _=ax.imshow(img, cmap=plt.get_cmap('gray'))
    _=ax.set_title(title)
    _=ax.axis('off');


image0, image1 = vortex()

fig, ax = plt.subplots(1,2, figsize=(16,8))
plot_grayscale(image0, ax=ax[0], title='image0')
plot_grayscale(image1, ax=ax[1], title='image1')

Downloading file 'data/pivchallenge-B-B001_1.tif' from 'https://gitlab.com/scikit-image/data/-/raw/2cdc5ce89b334d28f06a58c9f0ca21aa6992a5ba/pivchallenge/B/B001_1.tif' to '/Users/Bethun/Library/Caches/scikit-image/0.23.2'.
Downloading file 'data/pivchallenge-B-B001_2.tif' from 'https://gitlab.com/scikit-image/data/-/raw/2cdc5ce89b334d28f06a58c9f0ca21aa6992a5ba/pivchallenge/B/B001_2.tif' to '/Users/Bethun/Library/Caches/scikit-image/0.23.2'.
No description has been provided for this image

Exercise 12-1: The iLK algorithm uses no smoothness constraint for regularization, imposing only the constant flow constraint. In a sentence of two, explain why for rough (random field) image the ikL algorithm with no smoothness constraint is a good choice.

Answer: The iLK algorithm without smoothness constraint is advantageous for rough, random field images because it allows for the preservation of fine details and intricate structures without enforcing artificial smoothness, thus better capturing the inherent irregularities and complexity in images.

You will now compute the optical flow components using iLK algorithm provided by skimage.registration.optical_flow_ilk. The code in the cell below does the follow operations:

  1. Compute the horizontal and vertical optical flow between the two images.
  2. Compute the Euclidean norm of the flow field.
  3. Plot the results.

Execute this code and examine the results.

In [3]:
# Compute the optical flow
v, u = optical_flow_ilk(image0, image1, radius=15)
# Compute flow magnitude
norm = np.sqrt(u ** 2 + v ** 2)

fig, ax = plt.subplots(2,2, figsize=(16,16))
plot_grayscale(v, ax=ax[0,0], title='Vertical flow')
plot_grayscale(u, ax=ax[0,1], title='Horizontal flow')
plot_grayscale(norm, ax=ax[1,0], title='Norm of flow')
ax[1,1].axis('off')
Out[3]:
(0.0, 1.0, 0.0, 1.0)
No description has been provided for this image

With the flow field computed, it is instructive to display the examine the vector flow field. The code in the cell below does just this using the quiver plot function matplotlib.pyplot.quiver. Execute this code and examine the result.

In [4]:
# Quiver plot arguments
nvec = 20  # Number of vectors to be displayed along each image dimension
nl, nc = image0.shape
step = max(nl//nvec, nc//nvec)

y, x = np.mgrid[:nl:step, :nc:step]
u_ = u[::step, ::step]
v_ = v[::step, ::step]

plt.imshow(norm)
plt.quiver(x, y, u_, v_, color='r', units='dots',
           angles='xy', scale_units='xy', lw=3)
plt.title("Optical flow magnitude and vector field")
plt.axis('off')

plt.show()
No description has been provided for this image

Exercise 12-2: Examine the plots of the horizontal flow, vertical flow, Euclidean norm of flow and the quiver plot of the flow field. Keep in mind that for the gray-scale images of flow a bright (light) value is strongly positive and a dark value is strongly negative. In a few sentences, answer the following questions.

  1. How does the flow field change as you move away from the axis of rotation?
  2. Are the different figures showing the flow field consistent and why? Hint; compare the vector flow to the images of horizontal and vertical flow as well as the norm of the flow.

Answers:

  1. We observe that as we move away from the axis of rotation in the flow field, the magnitude of the flow tends to increase, indicating higher velocities. Near the axis of rotation, the flow vectors are relatively uniform and aligned with the axis, while farther away, they exhibit more divergence, representing the rotational motion spreading across the image. This pattern is visible in both the magnitude and direction of the flow field, as well as in the quiver plot.
  2. Yes, the different figures showing the flow field are consistent. The quiver plot of the flow field displays vectors representing the horizontal and vertical components of the flow. These components are also individually visualized in the plots of horizontal and vertical flow images. Additionally, the norm of the flow image combines the magnitudes of the horizontal and vertical flow components into a single image, which corresponds to the overall intensity of the flow represented by the arrows in the quiver plot.

TV-L1 Algorithm¶

You will now work with an improved version of the TV-L1 algorithm. Here you will examine the optical flow between two time-sequential images. This image data set comes with ground truth.

Execute the code in the cell below to load the

In [5]:
# Load the sequence of images and displacement
image0, image1, disp = stereo_motorcycle()

The code in the cell below displays the two images, along with the Euclidean norm of ground truth displacement. Execute this code and examine the results.

In [6]:
image0 = rgb2gray(image0)
image1 = rgb2gray(image1)

fig, ax = plt.subplots(2,2, figsize=(16,12))
plot_grayscale(image0, ax=ax[0,0], title='image0')
plot_grayscale(image1, ax=ax[0,1], title='image1')
plot_grayscale(disp, ax=ax[1,0], title='Displacement ground truth')  
ax[1,1].axis('off')
Out[6]:
(0.0, 1.0, 0.0, 1.0)
No description has been provided for this image

The lighter the pixel in the displacement image the greater the displacement.

Another way to view the displacement between these images is by applying an optical trick. A 3-channel (RGB) image is created from the gray-scale images. The second (displaced) image is placed in the first channel. The first (reference) image is placed in the second and third channel. In this way, displacements can be seen as red-green shadows. Execute the code in the cell below and examine the result.

In [7]:
def DisplacementTORGB(image0, image1):
    nr,nc = image0.shape
    seq_im = np.zeros((nr, nc, 3))
    seq_im[..., 0] = image1
    seq_im[..., 1] = image0
    seq_im[..., 2] = image0
    return(seq_im)

plot_grayscale(DisplacementTORGB(image0, image1))
No description has been provided for this image

The code in the cell below performs the following operations:

  1. Computes the optical flow between the two images using the TV-L1 algorithm provided by the skimage.registration.optical_flow_tvl1 function.
  2. Computes norm of the flow field vectors.
  3. Computes the pixel-by-pixel end-point error, filter for infinite values (from zero divides), and print the average end-point error.
  4. Display the vertical flow field, horizontal flow field, Euclidean norm of the flow, displacement ground truth, and end-point error field.

Execute the code and examine the results.

In [8]:
v, u = optical_flow_tvl1(image0, image1)

flow_norm = np.sqrt(u ** 2 + v ** 2)

endpoint_error = np.sqrt(np.square(flow_norm - disp))
endpoint_error[endpoint_error==np.Inf] = 0.0
print('Average end-point error = ' + str(np.mean(endpoint_error)))

fig, ax = plt.subplots(3,2, figsize=(16,18))
plot_grayscale(v, ax=ax[0,0], title='Vertical flow')
plot_grayscale(u, ax=ax[0,1], title='Horizontal flow')
plot_grayscale(flow_norm, ax=ax[1,0], title='Norm of flow')
plot_grayscale(disp, ax=ax[1,1], title='Displacement ground truth')  
plot_grayscale(endpoint_error, ax=ax[2,0], title='End-point error field')  
ax[2,1].axis('off')
ax[2,1].axis('off')

Average end-point error = 6.4592657
Out[8]:
(0.0, 1.0, 0.0, 1.0)
No description has been provided for this image

Exercise 12-3: Examine these results and answer the following questions in one or two sentences.

  1. This the estimated flow more vertical or horizontal, and how can you tell?
  2. What portion of the image exhibits the greatest error in the flow field and how can you tell.

Answers:

  1. The dominant flow appears more prominently in the horizontal direction (from above results), so the estimated flow is more horizontal.
  2. On inspecting the end-point error field, we observe that the error is higher in horizontal direction as the front and rear of the motorcycle appears brighter. This suggests that the front and rear of the motorcycle exhibits greatest error.

Next, execute the code in the code to display the original image first image, the warped second image, and the registered image created with the first image warped second image.

In [9]:
nr, nc = image0.shape
row_coords, col_coords = np.meshgrid(np.arange(nr), np.arange(nc), indexing='ij')
image1_warp = warp(image1, np.array([row_coords + v, col_coords + u]),mode='edge')

fig, ax = plt.subplots(2,2, figsize=(16,12))
plot_grayscale(DisplacementTORGB(image0, image0), ax=ax[0,0], title = 'Original image')
plot_grayscale(image1_warp, ax=ax[0,1], title = 'Warped image')
plot_grayscale(DisplacementTORGB(image0, image1_warp), ax=ax[1,0], title = 'Registered image')
ax[1,1].axis('off')
Out[9]:
(0.0, 1.0, 0.0, 1.0)
No description has been provided for this image

Exercise 12-4: Examine these results. Where are the largest registration errors and are these consistent with the flow errors you investigated in Exercise 12-3?

Answer: We observe that the front of the motorcycle in does not align correctly in the warped and registered images as compared to the original image. These misalignments indicate registration errors in the horizonal direction. By comparing the regions of high flow error with areas of misalignment in the registered image, we get the same conclusion as Exercise 12-3 that the estimated flow is more horizontal.

Next, you will repeat the steps of computing characteristics of the optical flow, but using the a larger parameter value, $\lamda$, 150, vs. the default value of 15 used in the first example. The larger parameter value favors brightness continuity. Execute the code and examine the results.

In [10]:
v, u = optical_flow_tvl1(image0, image1, attachment=150)

flow_norm = np.sqrt(u**2 + v**2)

endpoint_error = np.sqrt(np.square(flow_norm - disp))
endpoint_error[endpoint_error==np.Inf] = 0.0
print('Average end-point error = ' + str(np.mean(endpoint_error)))

fig, ax = plt.subplots(3,2, figsize=(16,18))
plot_grayscale(v, ax=ax[0,0], title='Vertical flow')
plot_grayscale(u, ax=ax[0,1], title='Horizontal flow')
plot_grayscale(flow_norm, ax=ax[1,0], title='Norm of flow')
plot_grayscale(disp, ax=ax[1,1], title='Displacement ground truth')  
plot_grayscale(endpoint_error, ax=ax[2,0], title='End-point error field')  
ax[2,1].axis('off')
ax[2,1].axis('off')

Average end-point error = 5.4036055
Out[10]:
(0.0, 1.0, 0.0, 1.0)
No description has been provided for this image

Exercise 12-5: Compare these results to the results with smaller value of $\lambda$ and answer the following these questions in one or two sentences.

  1. Comparing the average end-point errors, which value of $\lambda$ provides a better result?
  2. What differences can you see in the flow field images and is this expected given the value change in $\lambda$?
  3. What differences can you see in the end-point error fields between the two values of $\lambda$?

Answers:

  1. On comparing the average end-point errors, we observe that the smaller value of lambda = 15 provides a better result than when lambda is larger (150).
  2. We observe that with a larger lambda value (150), the images exhibit smoother and more globally consistent flow patterns, having overall brightness continuity. However, with a smaller lambda value (15), the images display more localized variations and details, as it captures finer structures in the flow.
  3. We observe that with a smaller lambda value (15), the end-point error field contains higher error values, especially in the front and rear of the motorcycle image. But, with a larger lambda value (150), the end-point error field shows lower error values overall, and also has greater emphasis on brightness continuity.

Finally, execute the code in the code in the cell below to display the registration of the first image with the warped second image.

In [11]:
row_coords, col_coords = np.meshgrid(np.arange(nr), np.arange(nc), indexing='ij')
image1_warp = warp(image1, np.array([row_coords + v, col_coords + u]),mode='edge')

fig, ax = plt.subplots(2,2, figsize=(16,12))
plot_grayscale(DisplacementTORGB(image0, image0), ax=ax[0,0], title = 'Original image')
plot_grayscale(image1_warp, ax=ax[0,1], title = 'Warped image')
plot_grayscale(DisplacementTORGB(image0, image1_warp), ax=ax[1,0], title = 'Registered image')
ax[1,1].axis('off')
Out[11]:
(0.0, 1.0, 0.0, 1.0)
No description has been provided for this image

These results are not particularly different from the first set. This is expected from the minimal change in average end-point error. While not great, you should be able to see that the error is in fact a bit less around the front wheel of the motorcycle.

Copyright 2022, 2024, Stephen F Elston. All rights reserved.¶

In [ ]: